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ABSTRACT 

We investigate the properties and evolution of accretion tori formed after the coalescence of 
two compact objects. At these extreme densities and temperatures, the accreting torus is cooled 
mainly by neutrino emission produced primarily by electron and positron capture on nucleons 
(/? reactions). We solve for the disc structure and its time evolution by introducing a detailed 
treatment of the equation of state which includes photodisintegration of helium, the condition of 
/3-equilibrium, and neutrino opacities. We self-consistently calculate the chemical equilibrium in 
the gas consisting of helium, free protons, neutrons and electron-positron pairs and compute the 
chemical potentials of the species, as well as the electron fraction throughout the disc. We find 
that, for sufficiently large accretion rates (M > IOMq/s), the inner regions of the disk become 
opaque and develop a viscous and thermal instability. The identification of this instability might 
be relevant for GRB observations. 

Subject headings: accretion, accretion discs - black hole physics - gamma rays: bursts - neutrinos 



1. Introduction 

Gamma-Ray Bursts are commonly thought to 
be produced in relativistic ejecta that dissipate en- 
ergy by internal shocks however alternative ideas 
based on the Poynting flux dominated jets are also 
being proposed (for a review see e.g. Piran 2005; 
Meszaros 2006; Zhang 2007). 

The enormous power released during the gamma- 
ray burst explosion indicates that a relativistic 
phenomenon must be involved in creating GRBs 
(Narayan et al. 1992). The merger of two neutron 
stars or of a neutron star and a black hole (e.g. NS- 
NS or NS-BH) has been invoked e.g. by Paczyhski 
(1986); Eichler et al. (1989); Paczyhski (1991); 
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Narayan, Paczyhski & Piran (1992), as well as the 
collapse of a massive star, the so called "collap- 
sar" scenario (e.g., Woosley 1993 and Paczyhski 
1998). In both cases a dense, hot accretion disk 
is likely to form around a newly born black hole 
(Witt et al. 1994). In the collapsar scenario, the 
collapsing envelope of the star accretes onto the 
newly formed black hole, while a transient debris 
disk is formed when the NS-NS or NS-BH binary 
merges (see e.g. Ruffert et al. 1997 for numerical 
simulations) . 

The durations of GRBs, which range from mil- 
liseconds to over a thousand of seconds, are dis- 
tributed in two distinct peaks defining two main 
GRB classes: short (< 2 sec) and long (> 2 sec) 
bursts (Kouvelietou et al. 1993). For some long 
bursts, signatures of an accompanying supernova 
explosion have been detected in the afterglow spec- 
tra (Stanek et al. 2003; Hjorth et al. 2003), which 
strongly favors the "collapsar" interpretation for 
their origin. Furthermore, the GRB positions in- 
ferred from the afterglow observations are consis- 
tent with the GRBs being associated with the star 
forming regions in their host galaxies. For short 
bursts, several pieces of evidence from the analy- 
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sis of the Swift satellite and follow-up observations 
(Gehrels et al. 2005; Fox et al. 2005; Villasenor et 
al. 2005) argue in favor of a binary merger model 
(Hjorth et al. 2005; Berger et al. 2005). 

In the merger scenario, the duration of the 
burst is comparable to the viscous timescale of the 
accretion disc whereas in the collapsar scenario, 
the external reservoir of stellar matter can feed 
the accretion torus for a much longer timescale. 
In general, accretion discs powering GRBs are ex- 
pected to have typical densities of the order of 
10 10-12 g cm~ 3 and temperatures of 10 11 K within 
10-20 Schwarzschild radii (R s = 2GM/c 2 ). Thus, 
accretion proceeds with rates of a fraction to sev- 
eral solar masses per second. In this "hyper- 
accreting" regime, photons become trapped and 
are not efficient at cooling the disc. Neutrinos, 
however, are produced by weak interactions in the 
dense and hot plasma, releasing the gravitational 
energy of the accretion flow. These discs go under 
the name of "neutrino-dominated accretion flows" , 
or NDAFs. 

Over the last several years a number of stud- 
ies have investigated the structure of these discs 
(Popham, Woosley & Fryer 1999; Narayan, Piran 
& Kumar 2001; Kohri & Mineshige 2002; Di Mat- 
teo, Perna & Narayan 2002; Surman & McLaugh- 
lin 2004; Kohri, Narayan & Piran 2005; Chen & 
Beloborodov 2006; Gu, Liu & Lu 2006; Liu, Gu 
& Lu 2007). These models have employed the 
customary approximation of one-dimensional hy- 
drodynamics (Shakura & Sunyaev 1973), where 
the effects of MHD viscous stresses are described 
by the dimensionless parameter a, but have been 
limited to the steady-state approximation of con- 
stant M. Such an assumption is a good approx- 
imation when considering the collapsar scenario, 
where the burst duration is much longer than the 
viscous time, due to the continuous replenishing of 
the disc by the collapsing star envelope. However, 
even in this scenario, recent observations suggest 
that that engines are "long-lived" (past the torus 
feeding phase), requiring a time-dependent com- 
putation. In the merger scenario, on the other 
hand, a time-dependent calculation is necessary 
even for modeling the prompt phase of the burst, 
since the duration is set by the viscous timescale 
of the disc. 

Recently, a fully time-dependent calculation of 
the structure of such accretion discs has been pre- 



sented by Janiuk et al. (2004). The model was 
suitable for the torus being a results of cither the 
gravitational collapse of a massive stellar core or 
the compact binary merger. However the struc- 
ture and evolution of the disk (like in many of the 
early calculations) was calculated under a num- 
ber of simplifying assumptions for the composition 
and the equation of state of the accreting matter. 
In this paper we improve upon our earlier results 
in several ways. Besides the requirement of time- 
dependent calculations, the high density and tem- 
perature regime in which the accreting gas lies, im- 
plies that both multi-dimensional numerical and 
semi-analytic calculations for such flows need to 
include the detailed microphysics. This includes 
photodisintcgration of nuclei, the establishment 
of statistical equilibrium, ncutronization, and the 
effects of neutrino opacities in the inner regions. 
Here, we introduce a detailed treatment of the 
equation of state, and calculate self-consistently 
the chemical equilibrium in the gas that consists 
of helium, free protons, neutrons and electron- 
positron pairs. We compute the chemical poten- 
tials of the species, as well as the electron fraction 
throughout the disk using the assumption of the 
equilibrium between the beta processes. Our EOS 
equations include, self-consistently, the contribu- 
tion of the neutrino trapping to the beta equilib- 
rium. Another important addition compared to 
our previous work (Janiuk et al. 2004) is the inclu- 
sion of photodisintegration of helium. The pres- 
ence of this term can affect the energy balance in 
the inner, opaque (to neutrinos as well as pho- 
tons) region of the flow and, as it will be shown, it 
eventually produces a thermal and viscous insta- 
bility in those regions. This is especially relevant 
since the GRB phenomenology requires a variable 
energy output. 

Other time-dependent disc studies of binary 
mergers or collapsars have been performed in 
2D using hydrodynamical simulations (e.g. Mac- 
Fadyen & Woosley 1999; Ruffert & Janka 1999; 
Lee et al. 2002; Rosswog et al. 2004; Lee, 
Ramirez-Ruiz & Page 2005) and, most recently, 
in 3-D simulations (Setiawan et al. 2005). Also, 
MHD simulations of the GRB central engine have 
been performed, showing that the magnetic field 
possibly plays an important role in the genera- 
tion of a GRB jet (Proga et al. 2003; Fujimoto 
et al. 2006). The advantage of our calculations 
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is that, whilst including all the relevant physics 
to calculate the equation of state, the structure 
and stability of the accretion disc, we are able to 
study a much larger range of parameter space and 
allow our calculations to evolve beyond what can 
be reached in higher dimensional calculation and 
comparable to at least the short-burst durations. 

The paper is organized as follows. In Sec- 
tion [3 we describe the basic assumptions of the 
model and the method used in the initial sta- 
tionary and subsequent time-dependent numerical 
simulations. In Section [3] we discuss the structure 
of the hyperaccreting disc for various values of the 
initial accretion rate, and study the time evolu- 
tion of its density and temperature, as well as the 
resulting neutrino lightcurve. We also discuss the 
physical origin of the instabilities in the disc, and 
we compare our model and results with the recent 
2D and 3D simulations. We summarize our results 
in Section [U 

2. Neutrino-cooled accretion disks 

In this Section, we describe how we improve 
upon our previous time-dependent calculation ( Ja- 
niuk et al. 2004) by computing self-consistently 
the equation of state of the extremely dense mat- 
ter by solving the balance of the (3 reaction rates. 
This allows us to determine the chemical poten- 
tials of electrons, protons and neutrons, as well as 
the electron fraction, in the initial disc configura- 
tion and throughout its evolution. 

2.1. Initial disc configuration: 1-D Hydro- 
dynamics 

We start by considering a steady-state model 
of an accretion disc around a Schwarzschild black 
hole - formed as a remnant structure either after 
a compact binary merger, or in a collapsar after 
the birth of a black hole (for a recent calculation 
in Kerr spacetime see Chen & Beloborodov 2006). 
Throughout our calculations we use the vertically 
integrated equations and hence derive a vertically 
averaged disc structure. We write the surface den- 
sity of the disk as £ = Hp, where p is the density 
and where the disk half thickness (or disk height) 
is given by H — c s /Qk- Here the sound speed 
is defined by c s = yP/p and Ok = \/GM/r 3 
is the Keplerian angular velocity with P the to- 
tal pressure. We note that, at very high accretion 



rates, the disc becomes moderately geometrically 
thick (H ~ 0.5r) in regions where neutrino cooling 
becomes inefficient and advection dominates. Our 
'slim disk' approximation neglects terms ~ (H/r) 2 
and assumes that the fluid is in Keplerian rota- 
tion. For the disc viscous stress we use the stan- 
dard a viscosity prescription of Shakura & Sun- 
yaev (1973) where the stress tensor is proportional 
to the pressure: 



-aP. 



(1) 



We adopt a value of a = 0.1 

We set the inner radius of the disc at 3 R$ , while 
the outer radius is at 50 Rs- The initial mass of 
such a disc is about O.35M0 for an accretion rate 
M = 1 Mq/s. Throughout the calculations we 
adopt a black hole mass of M = 3M Q . 

2.2. The equation of state 

We assume that the torus consists of helium, 
electron-positron pairs, free neutrons and protons. 
The total pressure is contributed by all particle 
species in the disc, and the fraction of each species 
is determined by self-consistently solving the bal- 
ance of the beta reaction rates. In the equation 
of state we take into account the pressure due to 
the free nuclei and pairs, helium, radiation and the 
trapped neutrinos: 



P = Pnucl + Pi 



He 



P* 



ad 



(2) 



The component P n uci includes free neutrons, pro- 
tons, and the electron-positron pair gas in beta 
equilibrium: 



P P , 
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with 



Pi = 
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where i*k are the Fermi-Dirac integrals of the or- 
der k, and rj , r] p and r] n are the reduced chemi- 
cal potentials of electrons, protons and neutrons 
in units of kT, respectively (where rji = pi/kT, 
also known as the degeneracy parameter, where pi 
the standard chemical potential) calculated from 
the chemical equilibrium condition (§ 12. 3[) . The 
reduced chemical potential of positrons is rj c+ = 



3 



—77c — 2/P e and the relativity parameters of the 
species i are defined as $ = kT/rmc 2 . 

Under the physical conditions in the torus, 
helium is generally non-relativistic and non- 
degenerate; therefore, its pressure is given by: 



■ffee = n ac kT, 



(5) 



where nn is the number density of helium. This 
is defined as: 



(6) 



and the fraction of free nucleons is given by 
X nuc = 295.5p- 3/4 T 1 9 1 /8 exp(-0.8209/T n ), (7) 

with Tu the temperature in unit of 10 11 K (e.g. 
Qian & Woosley 1996; Popham et al. 1999). 
The radiation pressure is given by: 
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When neutrinos become trapped in the disc, the 
neutrino pressure is non-zero. Following the treat- 
ment of photon transport under the two-stream 
approximation (Popham & Narayan 1995; Di Mat- 
teo et al. 2002), we have 
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where r s is the scattering optical depth due to 
the neutrino scattering on free neutrons and pro- 
tons and T a ,„ c and t^ v are the absorptive opti- 
cal depths for electron and muon neutrinos, re- 
spectively (see§2.3). The contribution from tau 
neutrinos is the same as that from muon neutri- 
nos. These optical depths and neutrino absorption 
processes (which are the reverse of the emission 
processes) are discussed in more detail in the Ap- 
pendix. 

In the disc we have to consider both the neu- 
trino transparent and opaque regions, as well as 
the transition between the two. In the trans- 
parent case, the neutrinos are not thermalizcd 
and the chemical potential of neutrinos is negli- 
gible. On the other hand, when neutrinos are 



totally trapped, the chemical equilibrium condi- 
tion yields: /i e + A*p = A*n + Hv The chemi- 
cal potential of neutrinos is a parameter depend- 
ing on how much neutrinos and anti-neutrinos are 
trapped, and assuming that the number densities 
of the trapped neutrinos and anti-neutrinos are 
the same, can be set to zero. In order to de- 
termine the distribution function of the partially 
trapped neutrinos, in principle one should solve 
the Boltzmann equation. To simplify this prob- 
lem, we use here a "gray body" model, and we 
introduce a blocking factor b = J2i= e a t ^ to de- 
scribe the extent to which neutrinos are trapped 
(see e.g. Sawyer 2003). In terms of this factor, we 
write the distribution function of neutrinos as 

&*(P) = — ( zbnXI = b ^> (0 < &i < !)■ 
exp(pc/Ki ) + 1 

(10) 

This simplified assumption is consistent with the 
two-stream approximation which we adopt here 
(Eq. E}. 

2.3. Composition and chemical equilib- 



The equilibrium state of the gas in the accret- 
ing torus is completely determined by the chem- 
ical potentials of neutrons, protons and electrons 
(r] n , 7] p , r) e ), and the trapping factor of neutrinos 
(6) which is related to the optical depths of neu- 
trinos (cf. Eq. 9). For a given baryon number 
density, rib, temperature T, and a value for accre- 
tion rate M and viscous constant a, the chemical 
potentials, or equivalently the ratio of free protons 
x = rip /rib, are determined from the condition of 
equilibrium between the transition reactions from 
neutrons to protons and from protons to neutrons. 
These reactions are: 



n + v c 

+ 



p + e 

p + D c — > n + e 
p + e~ + v e — > n 
n + e + — ► p + v c 
n —>■ p + e~ + v c 
n + v c — > p + e~ 



(11) 
(12) 

(13) 

(14) 

(15) 
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Therefore we have to calculate the ratio of protons 
that will satisfy the balance: 
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The reaction rates are the sum of forward and 
backward rates and are given in the Appendix (see 
also Kohri, Narayan and Piran 2005). 

These are supplemented by two additional con- 
ditions: the conservation of the baryon number, 
n n + n p = rib x X nuc , and charge neutrality (Yuan 
2005): 

p + ^c, (18) 



n P + 



which says that the net number of electrons is 
equal to the number of free protons plus the num- 
ber of protons in helium: 



(1 - X DUC )-^-. 



(19) 



The number density of fermions under arbitrary 
degeneracy is determined by the following equa- 
tions: 
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Finally, the electron fraction is defined as: 

_ n e - - n e + 



(21) 



(Note that this is different from Y c = 1/(1 + 
n n /n p ), which is only valid for free n-p-e gas.) 

2.4. Neutrino cooling 

The processes that are responsible for the neu- 
trino emission in the disc are electron-positron pair 
annihilation (e _ + e + — > V\ + V\), bremsstrahlung 
(n + n — > n + n + v\ + &{), plasmon decay (7 — > 
Vq. + v c ) and URCA process (reactions [TT1 [T4l and 
I15[) . The first two processes produce neutrinos of 
all flavors, while the other produce only electron 
neutrinos and anti-neutrinos. 

The cooling rate due to pair annihilation is ex- 
pressed as: 



<7c+c 



(22) 



where pio is the baryon density in units of 10 10 
g/cm 3 and T\\ is temperature in units of 10 11 K. 

The cooling rate due to the plasmon decay (in 
erg/cm 3 /s) is: 
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where 7p = 5.565 x 10" V^ 2 + 3(/x e /fcT) 2 )/3. 

The cooling rate due to the URCA reactions is 
given by the three emissivities: 

<?urca = Qp+e~ — >n+v e ' Q , n+e+— >p+P e t 5n— »p+e _ +S e ■ 

(25) 

The emissivities are given in the Appendix. 

Note that the blocking factor of the trapping 
neutrinos is used only for the emissivities of the 
URCA reactions. For simplicity, we neglect the 
blocking effects of neutrinos when calculating the 
emissivities for the electron-positron pair annihi- 
lation. Two reasons make this approximation rea- 
sonable: the emissivities for the electron-positron 
pair annihilation is much smaller than those of the 
URCA reactions, and the electron-positron pair 
annihilation does not change the electron fraction 
which sensitively affects the EOS. 

Each of the above neutrino emission process has 
a reverse process, which leads to neutrino absorp- 
tion. These are given by Equations 021 [13] and 
[TBI Therefore we introduce the absorptive optical 
depths for neutrinos given by: 



H 

o 



(26) 



where absorption of the electron neutrinos is de- 
termined by: 

asm ~r ybrems 

(27) 



and for the muon neutrinos: 
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(28) 



where the cooling rates for all three neutrino fla- 
vors are calculated by means of Fermi-Dirac inte- 
grals and are given in the Appendix. 

The cooling rate due to nucleon-nucleon bremsstrahlung T 
(in erg/cm 3 /s) is given by: 



In addition, the free escape of neutrinos from 
the disc is limited by scattering. The scattering 
optical depth is given by: 



?br, 
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where C s 



(4(Cy - l) 2 + 5a 2 )/24, C s , n = (1 + 



5a 2 )/24, CV = l/2 + 2sm 2 c , with a = 1.25 and 
sin 2 6c = 0.23. 

The neutrino cooling rate is then given by 

laT 4 
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(30) 
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2.5. Energy and momentum Conservation 

The hydrodynamic equations we solve to cal- 
culate the disc structure are the standard mass, 
energy and momentum conservation. 

Making use of the standard disk equations, the 
vertically integrated viscous heating rate (per unit 
area) over a half thickness H is given by: 



3GMM 
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(31) 



where the Newtonian boundary condition is as- 
sumed: f(r) = 1 — \Jr m i n /r. Note that in the 
time-dependent calculations, instead of Eq. [311 we 
will solve the viscous diffusion equation (Eq. [44]) . 

Using mass and momentum conservation M = 
AitpRHvr w QnvpH where v r w (3f)/(2r) and 
v = (2Pa)/(3/jf2) is the kinematic viscosity. The 
viscous heating rate can be written in terms of a, 

3 



-aVLHP. 



(32) 



Cooling in the disc is due to advection, radia- 
tion and neutrino emission. The advective cooling 
in a stationary disc is determined approximately 
as: 



aPHT 



(33) 



where q a( jv oc d In S/d In r oc (d In T/d In r — (T^ — 
l)dhi p/d\nr) m const and we adopt the value 
of 1.0. The entropy density S is the sum of four 
components: 



S — <5nucl ~F Shc F <SVad + S v 



(34) 



The entropy density of the gas of free protons, 
neutrons and electron-positron pairs is given by: 



5*nucl — S e - + S c + F S p F S n 



where 
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and 
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is the energy density of electrons, positrons, pro- 
tons or neutrons, Pi is their pressure, given by 
Equation |4l m are the number densities and rji are 
the chemical potentials. 

The entropy density of helium is given by: 



Sac 



"-Ho 
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for 7lHc > 0. 

The entropy density of radiation is 



s 



rad 



I ^rad 



while for neutrinos we have 

p 

S = 4— 
kT' 



(38) 



(39) 



(40) 



In the initial disc configuration we assume that 
<Zadv is approximately constant and of order of 
unity, but in the subsequent time-dependent evo- 
lution the advection term is calculated with the 
appropriate radial derivatives. 

For the case of photon and electron-positron 
pairs in the plasma the radiative cooling is equal 
to: 

3P rad c llaT 4 
= — = (41) 
where we adopt the Rosseland-mean opacity k = 
0.4 + 0.64 x 10 23 pT- 3 [cm 2 g -1 ]. 

An important term in the cooling and heating 
balance in the disc is due to photodisintegration 
of a particles, with rate: 



^5photo — Qphoto-^ 



where 



<? P hoto = 6.28 x 10 28 p 10 v r 



dX n 
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(43) 



and X nuc is given by Equation [7] Finally, in order 
to calculate the initial stationary configuration, we 



solve the energy balance: F t( 
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2.6. Time evolution 

After solving for the initial disc configuration, 
we allow the density and temperature to vary with 
time. We solve the time-dependent equations of 
mass and angular momentum conservation in the 
disc: 



3r 



SE _ ia_ 

dt r dr 
and the energy equation 
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where x = (P ~ Pvadi/P- The cooling term Q~ 
consists of radiative and neutrino cooling, given by 
Equations (|4"Tj) and ([50]) . Advection is included in 
the energy equation via the radial derivatives. The 
cooling term due to photodisintegration of helium 
now must be proportional to the full time deriva- 
tive of X nuc (cf. Eq. 



Qphoto OC V, 



2.7. Numerical method 



dX n 



Of 



(46) 



The initial configuration of the disc is calcu- 
lated by means of the Newton-Raphson method, 
iterated with the hydrostatic equilibrium condi- 
tion. We interpolate over the matrix of pre- 
calculated results for the equation of state (pres- 
sure and entropy) and neutrino cooling rate ( 
the number of points is 1024x1024). The Fermi- 
Dirac integrals are calculated using the mixture of 
Gauss-Legendre and Gauss-Laguerre quadratures 
(Aparicio 1998). 

Having determined the initial radial profiles of 
density and temperature, as well as the other 
quantities at time t = 0, we start the time evolu- 
tion of the disc. We solve the set of Equations (j44[) , 
(|46| and (|46]) using the convenient change of vari- 
ables y — 2r x / 2 and 5 = ?/E, at fixed radial grid, 
equally spaced in y (see Janiuk et al. 2002 and ref- 
erences therein) . The number of radial zones is set 
to 200, which we found to be an adequate resolu- 
tion. After determining the solutions for the first 
100 time steps by the fourth-order Runge-Kutta 
method, we use the Adams-Moulton predictor- 
corrector method with an adaptive time step. The 



code used an explicit communication model that 
is implemented with the standardized MPI com- 
munication interface and can be run on multipro- 
cessor machines. 

We choose the no-torque inner boundary con- 
dition, Ei n = T; n = (see Abramowicz & Kato 
1989). The outer boundary of the disc is done by 
adding an extra "dead-zone" to the computational 
domain, which accounts for the disk expansion and 
conservation of angular momentum. 

3. Results 

We first analyze the pressure, entropy and neu- 
trino cooling rate distributions for a given tem- 
perature and baryon density in the gas. Then, 
we show the disc structure for a converged static 
disc model and finally we show examples of time 
evolution of the neutrino luminosity, density, tem- 
perature and electron fraction for given sets of pa- 
rameters. 

3.1. EOS solutions for a given tempera- 
ture and density 

In Figure Q] and [2] we plot the results of the nu- 
merically calculated equation of state for the hot 
and dense matter. The plots show the dependence 
of the electron fraction, pressure, entropy and neu- 
trino cooling rate on temperature and density, re- 
spectively. 

In the upper panels, we show the neutrino cool- 
ing rate. At low temperatures, below T — m c c 2 ~ 
5 x 10 9 K, there are almost no positrons and free 
nucleons. Therefore the neutrino emission pro- 
cesses switch off, and the cooling of the gas is 
either due to advection, or, when the matter be- 
comes transparent to photons, radiative cooling 
overtakes. For larger temperatures, the neutrino 
emission rate increases up to the temperature of 
about ~ 5 x 10 11 K. For very high tempera- 
tures, the optical depths for neutrinos increase 
very rapidly (r cx T 5 , see Eq. (7) in Di Matteo 
et al. 2002). Therefore the neutrino cooling rate 
decreases at high temperatures (Eq. [30]) . On the 
other hand, for a given temperature (e.g. T ~ 10 11 
K), the neutrino cooling rate does not sensitively 
depend on density. It varies by one order of magni- 
tude in the range of 10 8 < p < 10 14 g/cm 3 , where 
the optical depth is t ~ 100. The middle panels 
of Figures Q] and [21 show the entropy and pressure 
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as a function of temperature and density. At low 
temperatures, the entropy of gas is not important. 
The highly degenerate electrons do not give contri- 
bution to the entropy, while they are a dominant 
term in the pressure, which is therefore indepen- 
dent of temperature up to T ~ 5 x 10 10 K. When 
the temperature increases, helium becomes disin- 
tegrated into free nucleons at energy comparable 
to the binding energy of helium, and after that the 
radiation (including photons and electron-positron 
pairs) contributes mainly to the total entropy and 
pressure. Therefore, both these quantities rise 
with temperature. At high densities, the entropy 
is dominated by neutrons. Finally, the electron 
fraction is shown in the bottom panel of Figures 1 
and 2. At low temperatures, the electron fraction 
is equal to 0.5, it then decreases sharply as the 
helium nuclei become disintegrated. As the tem- 
perature further increases, positrons appear as the 
electrons become non-degenerate. The positron 
capture again increases the electron fraction (see 
Figfl). 

The electron fraction changes significantly as 
a function of density for T > 10 W K (in FigEl 
T = 10 n iT). At low densities, the torus consists 
of free neutrons and protons and Y e is close to 0.5 
(see also Eq. 7). As density increases, Y c decreases 
to satisfy the beta-equilibrium among the free n-p- 
e gas. Above some density (when the temperature 
is high enough, e.g. for Fig. [2j pu c ~ 10 13 g cm -3 ) 
helium starts forming. Therefore Y c has a kink 
and stars rising steeply, asymptotically approach- 
ing 0.5 as the torus consists of plenty of ionized 
helium and some electrons to keep charge neutral- 
ity. 

3.2. The steady-state disc structure 

In Figures [3] and |4] we show the profiles of den- 
sity and temperature in the stationary accretion 
disk model for three accretion rates: 1 M©/s, 
10 Mq/s and 12 Mq/s. In general, the temper- 
ature and density profiles both increase inward. 
However, for M = 12Af /s, a distinct branch of 
solutions is reached, which appears different than 
the so-called "NDAF" branch (see Kohri & Mi- 
neshige 2002). The density and temperature pro- 
files for this high accretion rate differ also from 
what was found in previous work (Di Matteo et 
al. 2002; Janiuk et al. 2004). Due to a more 
detailed equation of state, in which we allow for 
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Fig. 1. — The dependence of the electron frac- 
tion (bottom), pressure (middle bottom), en- 
tropy (middle upper) neutrino cooling rate (up- 
per panel) on temperature, for the constant den- 
sity p = 10 12 g/cm 3 . The accretion rate is M = 
1 Mq/s. The pressure and neutrino cooling are in 
cgs units and the entropy is in units of /cb cm~ 3 . 
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Fig. 2. — The dependence of the electron frac- 
tion (bottom), pressure (middle bottom), entropy 
(middle upper) and neutrino cooling rate (upper 
panel) on density, for the constant temperature 
T = 10 11 K. The accretion rate is M = 1 M Q /s. 
The pressure and neutrino cooling are in cgs units 
and the entropy is in units of /cbcto -3 . 
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a partial degeneracy of nucleons and electrons as 
well as neutrino trapping, our solutions reach den- 
sities as high as 10 12 g/cm 3 in the innermost radii 
of the disc. The temperature in this inner disc part 
is in the range 4 x 10 10 — 1.25 x 10 11 K, depending 
on the accretion rate. For the hottest disk model, 
a local peak in the density forms around 7 — 8Rs , 
while below that radius the density decreases. Be- 
tween ~ 3.5 and 7 Rs, the plasma becomes much 
hotter and less dense than outside of this region. 
This means that the macroscopic state of the sys- 
tem is different here due to an abrupt change in 
the heat capacity. In order to check what is the 
reason for this transition, we investigate the pres- 
sure distribution in the disk. 

The profile of the pressure is shown in Fig. [5] 
The dominant term in the total pressure is due 
to the nucleons, while the radiation pressure (in- 
cluding electron-positron pairs) is always several 
orders of magnitude smaller. The neutrino pres- 
sure is large in the inner disc, once it gets opti- 
cally thick to neutrinos (i.e. for M > IOA/q/s). 
A significant contribution to the pressure is due 
to helium at densities high enough for helium to 
form, albeit at temperatures low enough such that 
its nuclei are not fully disintegrated. For the 
largest accretion rate shown, in the region of the 
temperature excess and inverse density gradient 
(3.5 — 7Rs), the total pressure distribution flat- 
tens. The helium pressure is now vanishingly small 
due to the complete photodisintegration, and the 
nuclear pressure is slightly decreased due to the 
composition change: smaller number density of 
neutrons and larger number density of protons. 
The substantial contribution to the pressure is now 
given by the neutrinos (large optical depths; see 
below) and radiation pressure (increased number 
of electron-positron pairs). From the comparison 
of Figures 02 [4] and the bottom panel of Figure 
03 it can be seen that the total pressure becomes 
locally correlated with temperature and anticor- 
related with density, thus consituting an unstable 
phase. 

In FigureOo]we show the neutrino optical depths 
due to scattering and absorption. The total opti- 
cal depth in the outer disc is typically dominated 
by scattering processes, while in the inner disc ab- 
sorption processes take over for very high accre- 
tion rates. For M = lM Q /s only the very in- 
ner disk radii have optical depth close to 1. For 
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Fig. 3. — The baryon density as a function of 
the disc radius, calculated in the stationary so- 
lution. The accretion rate is M = 1 M@/s (solid 
line), M = 10 M©/s (long dashed line) and M = 
12 Mq/s (short dashed line) . 
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Fig. 4. — The temperature as a function of the disc 
radius, calculated in the stationary solution. The 
accretion rate is M = 1 M Q /s (solid line), M 
1 Mq/s (long dashed line) and M = 12 M Q /s 
(short dashed line) . 
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Fig. 5. — The pressure components as a func- 
tion of the disc radius, calculated in the station- 
ary solution for the three accretion rate values: 
M = 1 M Q /s (upper panel), M = 10 M /s (mid- 
dle panel) and M = 12 M /s (bottom panel) 
The total pressure is marked by the solid 
line, and its components are: nuclear (gas) pres- 
sure (long dashed line), radiation pressure (short 
dashed line), helium pressure (dotted line) and 
neutrino pressure (dot-dashed line). 
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Fig. 6. — The neutrino optical depths due to scat- 
tering (t s , solid line) and absorption (r aje for elec- 
tron neutrinos, long dashed line, and T a .^ for muon 
neutrinos, short dashed line) as a function of ra- 
dius for M = IMq/s (upper panel), M = 10M o /s 
(middle panel) and for M = 12M Q /s (bottom 
panel). The sum of the three quantities is the 
total optical depth (r tot , dotted line). 



M = 12Af /s, in the radial strip of - 3.5 - 7R S 
the disk is optically thick with absorptive optical 
depth for electron neutrinos exceeding the scatter- 
ing term and reaching values of the order of 100. 

3.2.1. Composition and Chemical potentials 

In Figure [7] we show the distribution of the 
reduced chemical potentials of protons, electrons 
and neutrons throughout the disc. Reduced elec- 
tron chemical potentials much larger than unity 
(indicating strong electron degeneracy) are found 
in the inner disc parts for M = IOMq/s and 
M = 12M /s, whereas for 1 Mq/s electrons are 
only slightly degenerate. For the highest accretion 
rate, the maximum degeneracies correspond to the 
radius of the local peak in the density (cf. Fig. [3|) 
and the excess of helium number density (cf. Fig. 
EJ. Below this radius, the species become non- 
degenerate again, contributing to the increase of 
the electron fraction (cf . Fig. [9]) . In Figure [8] we 
plot the mass fraction of free nucleons as a function 
of radius for M = 12M /s, 10 and lM /s. As the 
Figure shows, in the outer regions, X nuc increases 
as the radius decreases, while the temperature and 
density increase (Fig. [3] and Fig. __]). Consistent 
with the behavior of Y e (Figj2]), X nuc subsequently 
turns around (decreases) at radii where the den- 
sity is high enough for significant helium forma- 
tion. This trend is reversed sharply for highest 
accretion rates, when the temperature in the disk 
is high enough (Fig_4j) for helium to be fully dis- 
sociated. In consequence, the number density of 
alpha particles increases at ~ 7— 12i?g and sharply 
decreases at lower radii. A similar, but far less pro- 
nounced fluctuation in X nuc is seen at smaller radii 
for the case of M = lOM /s. For smallest accre- 
tion rate, lM /s, there is no helium throughout 
the disk. 

In Figure _5] we show the radial distribution of 
the electron fraction throughout the disc for 1, 
10 and 12 M /s. For the case of 1 M /s (solid 
line), the electron fraction decreases inward in the 
disc as the electrons are captured by protons (in 
neutronization reactions). Once the electrons be- 
come non-degenerate, positrons appear, and the 
positron capture by neutrons again increases the 
electron fraction. For the hotter plasma (accre- 
tion rate of 10 Mq/s, dashed line), consistently 
with the behavior discussed for X nuc , helium nu- 
clei form as the density becomes high enough be- 
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Fig. 7. — The chemical potentials of neutrons 
(i] n , solid line), electrons (r]e, dashed line) and 
protons (?7p, dotted line) as a function of the 
disc radius, calculated in the stationary solution. 
The accretion rate is M — 1 M Q /s (triangles), 
M = 10 Mq/s (squares) and M = 12 Af /s (cir- 
cles). 
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Fig. 8. — The mass fraction of free nucleons 
as a function of radius for M = 1Mq/s (solid 
line), M = 10M Q /s (long dashed line) and M = 
12 Mq/s (short dashed line). 



low ~ 20i?s and Y e increases. For the accretion 
rate of 12 Mq/s there is the sharp decrease in Y e , 
at ~ 7 — 8R$, due to the sudden dissociation of he- 
lium. As helium is fully photo-dissociated, there is 
an almost equal number of neutrons and protons 
due to the balance of the electron and positron 
capture. This implies an electron fraction of 0.5. 
At the innermost radius, the temperature and den- 
sity drop due to the boundary condition, which 
affects the behaviour of both Y e and A nuc . 

3.2.2. Cooling and heating rates 

In Figure [TU] we plot the rates of viscous heat- 
ing, advection and cooling due to neutrino emis- 
sion and photo-dissociation in the stationary disc. 
The accretion rates are M = 1 Mq/s (upper 
panel), M = 10 Mq/s (middle panel), and M = 
12 Mq/s (lower panel). For the highest accretion 
rates, in the innermost disc the neutrino cooling 
rate decreases substantially with respect to the 
cooling by photodissociacion. This is because the 
neutrinos are trapped in the disc due to a large 
opacity The smaller the accretion rate, the less 
important is the neutrino trapping effect. This 
implies that for an accretion rate of < 10 Mq/s 
neutrinos can escape from the innermost disc. 

The advective term is a couple orders of magni- 
tude smaller than the other terms. The photodis- 
sociacion term is negligible for an accretion rate of 
1 Mq/s, since there is no helium in the whole disc, 
and Qphoto is equal to zero by definition. For an 
accretion rate of 10 Mq/s there is very little he- 
lium down to about 15-20 i?s, and therefore Q p hoto 
is much smaller than other terms. For the accre- 
tion rate of 12 Mq/s , down to 6 — 10i?s m the 
region of the disc of high density and maximum 
degeneracy, helium nuclei form. The nucleosyn- 
thesis of alpha particles leads to the plasma heat- 
ing instead of cooling, and therefore the relevant 
term in the energy balance has a negative value. 
Outward, above ~ 10i?s, there is some fraction 
of helium which can be photo-dissociated, so the 
cooling term due to this reaction is also important 
in the total energy balance. In the inner region 
helium is fully dissociated and Q p hoto is equal to 
zero, increasing again only near the inner bound- 
ary due to the local density increase and decrease 
of temperature. 
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Fig. 9. — The electron fraction as a function 
of the disc radius, calculated in the stationary 
solution. The accretion rate is M = 1 Mq/s 
(solid line), M = 10 M©/s (long dashed line) and 
M = 12 Mq/s (short dashed line). 
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Fig. 10. — The heating and cooling rates due 
to photodissociation and neutrino emission (solid 
lines) as a function of radius for M = lM /s (up- 
per panel), M = IQMq/s (middle panel) and for 
M = I2M0/S (bottom panel). The other terms 
are: cooling rate due to advection (long dashed 
line) and viscous heating rate (short dashed line). 



3.3. Stability analysis: instabilities at 
high- M 

The disc is thermally unstable if d log Q + jd log T > 
d\ogQ~ /dlogT. Then any small increase (de- 
crease) in temperature leads to a heating rate 
which is more (less) than the cooling rate, and 
as a consequence a further increase (decrease) of 
the temperature. The viscous instability, which 
appears when |>|f| q+_q- < 0, manifests itself 
in a faster (slower) evolution of an underdense 
(overdense) region. The instabilities can be con- 
veniently located in the surface density - temper- 
ature diagrams, in which the branch of thermal 
equilibrium solutions with a negative slope is not 
only unstable to the perturbations in the surface 
density, but it is also thermally unstable. 

In Figure [Tl] we show such stability curves for 
several radii in the disc. The criterion for a vis- 
cously stable disc is generally satisfied through- 
out the whole disc for M < 10M Q /s. How- 
ever, for larger accretion rates, there are unstable 
branches at the smallest radii. For M = 10M & /s, 
the disc becomes unstable below 5 i?s, while for 
M = 12M Q /s the instability strip is up to ~ 7Rs- 
Here helium is almost completely photodisinte- 
grated while the electrons and protons become 
non-degenerate again. For this high accretion rate, 
the electron fraction rises inward in the disc. Un- 
der these conditions, the energy balance is affected 
leading to the thermal and viscous instability, as 
demonstrated by the stability curves. This insta- 
bility will be discussed in more detail in Section 

S3J 

3.4. Time dependent solutions 

In this Section, we discuss how the tempera- 
ture, density, electron fraction and disk luminosity 
evolve with time. In Figures[T2land[T3lwe show the 
time evolution of density and temperature, when 
the initial accretion rate is 1 Mq/s. These quan- 
tities exponentially decrease with time: 

p = p (r) exp(-at) , (47) 

and 

T = T (r)exp(-6t) (48) 

where a ~ 1.9 and b w 0.085. The normalization 
of these relations depends on the radius, and for 
example for r = 6i?s it is po = 2.2 x 10 11 and 
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To = 3.5 x 10 10 . The exponential behaviour arises 
from the nature of energy equation (45). 

In Figure [14] we show the electron fraction as 
a function of time for several exemplary radial lo- 
cations in the disc, for the disc evolving from a 
starting accretion rate of 1Mq/s. The fraction Y e 
is smaller in the inner disc radii, while outward, 
the electron fraction is over half an order of mag- 
nitude higher. Altogether, during the evolution 
of the system, the electron fraction constantly in- 
creases with time throughout the disc. 

The time-dependent neutrino luminosity of the 
disc is given by: 

L v (t)= [ Q-(t)2nrdr (49) 

where Q~ is given by Equation (f30|) . 

In Figure [T5] we show an example of such a 
lightcurve, for our standard model parameters 
(M = 3M Q , a = 0.1), and i? max = 50i? S - The 
starting accretion rate was M s tart = 1 M /s. At 
this accretion rate neutrinos can already escape 
from the accretion disc at the beginning of the 
evolution. For higher initial accretion rates, e.g 
10 — V2Mq/s, neutrinos are trapped in the inner- 
most disc, and, as a consequence, the neutrino 
luminosity is lower at the initial stages of disc 
evolution until the accretion rate drops to about 
~ IMq/s. This result is qualitatively similar, al- 
beit it differs quantitatively, from what was ob- 
tained in Di Matteo et al. (2002) and Janiuk et 
al. (2004): in those calculations neutrino trap- 
ping was far more substantial even for a 'mod- 
erate' accretion rate of M > lM /s. The differ- 
ence arises from the fact that here we calculate the 
neutrino opacities using the f3 reaction efficiencies, 
self-consistently with the equation of state. 

For an accretion rate of lM Q /s, the solution 
does not reach the viscously unstable branch. Ini- 
tially, the disc contains almost no alpha particles 
(cf. Fig. [8|), which appear later on during the evo- 
lution and cooling of the plasma. The dynamical 
balance between the photodisintegration of helium 
and nucleosynthesis leads to an additional non- 
zero cooling/heating term in the energy equation 
and to only small amplitude flickering at the early 
stages of time-evolution. 

The situation is much more dramatic when the 
starting accretion rate is 12M Q /s. In this case a 
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Fig. 11. — The stability curves on the accretion 
rate vs. surface density plane, for several chosen 
radii in the disc: 3.39i?s (solid line), 3.81i?s (dot- 
ted line), 4.25i?s (short dashed line), 5.19i?s (long 
dashed line) and 8.60i?s (dot-dashed line). 
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Fig. 12. — The density as a function of time, for 
several chosen disc radii: 4.01, 6.04, 10.13, 20.7, 
35.02, and 45.19 Rs- The initial accretion rate is 
M = 1 Mq/s. 
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Fig. 13. — The temperature as a function of time, 
for several chosen disc radii: 4.01, 6.04, 10.13, 
20.7, 35.02, and 45.19 R s . The initial accretion 
rate is M = 1 M /s. 
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Fig. 14. — The electron fraction as a function 
of time, for several chosen disc radii: 4.01, 6.04, 
10.13, 20.7, 35.02, and 45.19 R s . The initial ac- 
cretion rate is M = 1 Mq/s. 



large disc strip is viscously and thermally unstable 
and the most violent instability takes place around 
and below 7 - 12i? s . 

In Figure HH] we show the behavior of the lo- 
cal accretion rate in the unstable disc, at sev- 
eral chosen locations within the instability strip. 
Near ~ 12i?s, the accretion rate varies due to 
the large and rapidly changing photodisintegra- 
tion term (locally, it can become larger than the 
neutrino cooling rate). 

This radius corresponds to the largest local 
value of the density of helium (cf. Figure O show- 
ing its starting model distribution), which is then 
being photodissociated. The photodissociacion 
process is the cause of the local rapid accretion 
rate changes. Then, inside from this highly vari- 
able strip, the accretion rate grows too fast to pre- 
serve the disc structure. This kind of behaviour 
occurs in the locally hotter and less dense region 
visible in the starting configuration e.g. in Figures 
[3]and|4l between 3.5 and 7 Rs- In this region the 
helium is already totally photodissociated. Due 
to the growing accretion rate all the material is 
rapidly accreted onto the black hole and the in- 
nermost strip of the disc empties. 

After the inner strip is destroyed, the outer 
parts can still accrete onto the center. As they 
approach the black hole, their temperature and 
density grow and the above situation can repeat 
several times, until the whole disc is completely 
broken into rings and destroyed. These later in- 
jections of energy, with timescales dictated by the 
viscous timescale of each ring, can produce en- 
ergy flares following the main GRB activity. Our 
results therefore provide another physical mech- 
anising for the flare model recently proposed by 
Perna et al. (2006). 

In Figure fT71 we show the neutrino lightcurve of 
the unstable disc. The instabilities due to photo- 
disintegration are reflected in oscillations of vari- 
able amplitude and millisecond timescale. This is 
of a particular interest if the neutrino annihila- 
tion provides the energy input for GRBs, however 
it should be pointed out that the oscillations ap- 
pearing in the presented lightcurve have a much 



1 ln addition to the gravitational instability in the outer parts 
of the disk, which was hinted by the calculations of Di 
Matteo et al. (2002) and confirmed by those of Chen & 
Beloborodov (2006). 
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smaller amplitude than the observed gamma ray 
variability. 
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Fig. 15. — The neutrino lightcurve, integrated 
over the disc surface. The initial accretion rate 
is M = 1 Mq/s. 
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Fig. 16. — The local accretion rate, as a function of 
time, for several chosen radial locations in the disc: 
6.86, 7.73 and 12.85 Rs- The starting accretion 
rate is M = 12 Mq/s. 



4. Discussion 

4.1. The unstable neutrino-opaque disc 

In our calculations we have shown that, for 
large accretion rates, the accreting torus becomes 
viscously and thermally unstable. We now discuss 
the physical origin of the instability. 

The unstable branch appears both in the 
steady-state solutions and in the subsequent time- 
dependent evolutions. In the steady-state case, 
for a chosen value of a constant accretion rate, 
this can be seen for instance by plotting the ra- 
dial profiles of density and temperature (cf. Figs. 
[3] and [H where the distinct branch is found for 
the innermost radii), as well as by looking at the 
stability curves for a range of accretion rates at 
a chosen disc radius (cf. Fig. QTJ where the un- 
stable inner disc radii exhibit a negative slope in 
the curve). In the time-dependent simulations, 
the unstable behavior is manifested by the highly 
variable accretion rate in certain strips of the disk 
and by the subsequent breaking of the disk inside 
from these variable strips (cf. Fig. IT5|) . The in- 
stability arises from the fact that the accretion 
rate rises locally too fast to prevent the disc strip 
from emptying, as the material is supplied from 
outer strips at much slower rate than it is accreted 
inwards. The disc evolves unstably on a viscous 
timescale, r v i S c = l/(Qif2) x (r/H) 2 ; for the radii 
shown in Figure [16j it is r visc w 0.05 s (note that 
the disc is rather thick, r/H ~ 2.5, and therefore 
the viscous and thermal timescales are close to 
each other). Theoretically, in order to find again 
a stable solution, the disc would have to increase 
the local accretion rate up to about several tens 
of Mq/s within one viscous timescale. However, 
this may not be possible if there is not enough 
material in the system to support much higher 
accretion rates during such violent oscillations. 
Therefore the system is unable to be stabilized 
and gets broken after a fraction of r v j sc . In addi- 
tion, the dynamical instability is the source of the 
flickering of the local accretion rate at the edge of 
the unstable strip. 

Let us now discuss in more detail the physical 
reason driving this instability. In the inner part of 
the disc (below r ~ 10i? s for M = 12M /s) there 
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Fig. 17. — The neutrino luminosity, as a func- 
tion of time The starting accretion rate is M = 
12 M Q /s. 
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Fig. 18. — The total pressure, as a function of 
time, for the chosen radial locations in the disc: 
7.73 and 6.86 i?s- The starting accretion rate is 
M = 12 M Q /s. 



are two important processes, both of which are 
incorporated in our equation of state: photodisin- 
tegration of helium and neutrino trapping. As it 
was already mentioned in Sec. 13.2.11 and can be 
seen from Fig. [51 below ~ 7 — 8i?s helium in this 
disc is completely photodisintegrated. This part 
of the disc is also opaque to neutrinos, as we show 
in Fig. [6] 

These two mechanisms competitively influence 
the electron fraction in the disc (cf. Sec. 13.21 Fig. 
[6] and Fig. [9]). Well outside the unstable strip, 
above ~ 20i?s, the electron fraction smoothly de- 
creases inwards as positrons appear because of the 
neutronization process. Then the scattering opti- 
cal depth for neutrinos becomes r s > 1, and the 
electron fraction increases again. After photodis- 
integration, the electron fraction decreases signif- 
icantly from almost 0.3 to much less than 0.1 due 
to electron capture. But again, when the disc 
becomes optically thick to absorption of electron 
neutrinos, the electron fraction gets higher and ap- 
proaches almost 0.5. 

The total pressure of sub-nuclear matter (cf. 
Fig [5|) is mainly contributed by electrons, and 
therefore it is influenced by the changes in the 
electron fraction. In the narrow range of radii 
(6.8 < r/Rs < 7.8), the pressure decreases due 
to photodisintegration. The sudden decrease of 
the pressure might drive the dynamical instabil- 
ity. (This picture is somewhat similar to that of 
the iron core collapse in the core collapse super- 
nova explosions: electron capture consumes most 
of the electrons and makes the EOS softer, conse- 
quently, it triggers the collapse of the iron core.) 
However, the transition from neutrino transpar- 
ent to opaque disc, and the increase of the elec- 
tron fraction due to the beta equilibrium (see also 
Yuan & Heyl 2005), are the reason for a steeper 
increase of the total pressure of the system. 

The same effect can also be observed in the 
time-dependent plot (Fig. [T5|) . in which we show 
the pressure changes in the characteristic radii 
of the unstable part of the disc (cf. Fig. [T6)) . 
The pressure decreases with time up to a radius 
R = 6.8i?S) since the temperature and density 
gradually drop, as well as the neutrino opacities, 
so the electron fraction gets smaller. Then, in a 
strip between ~ 6.8 and 12 i?s, the pressure rises 
with time: in fact, when a particles appear, the 
electron fraction rises and matter locally piles up, 
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thus increasing the pressure. 

At the border of these radii the disc breaks up, 
when the thermal-viscous instability induces an 
avalanche-like increase of the local accretion rate 
below ~ 8Rs- This happens because the increase 
in the pressure causes an excess in the local energy 
dissipation rate and the disc heats up, while at ad- 
jacent radii the pressure decreases and heating is 
insufficient. The system tries to compensate these 
temperature gradients by decreasing/ increasing 
the temperature in the outer/inner radius, respec- 
tively. But since in the unstable mode of the ther- 
mal balance this causes a further increase/decrease 
of density, the pressure drops further and the disc 
heats up in the outer radius, while cools down in 
the inner one. As long as it cannot find any sta- 
ble track of evolution, the emptying of the inner 
strip continues and finally the whole material is 
accreted towards the black hole or blown out. 

The radial extent of the unstable part depends 
on the initial accretion rate and in our model for 
M = 12M /s it is up to ~ 8R S , for M = 10M o /s 
it is up to ~ 5i?S: while for M = lM@/s it is below 
~ 3.5i?s- In the latter case, since the inner radius 
is located at ~ 3i?s, the instability hardly affects 
the disc. The extension of the instability strip de- 
pends also on the mass of the accreting compact 
object, and since for lower mass black holes the ac- 
cretion disc is generally denser, it reaches a density 
~ 10 12 g/cm 3 around ~ 15i?s- 

Above 25-30 Rs, where the plasma is already 
optically thin and the evolution is stable, both 
the pressure and the accretion rate smoothly drop 
with time. The dominant source of cooling of the 
disc in this region is the neutrino emission (ad- 
vective cooling decreases as the disc transits from 
neutrino opaque to transparent). The photodis- 
integration term (if non-zero), is usually by 1-2 
orders of magnitudes smaller than neutrino cool- 
ing, and in the inner disc, up to about 6.8 Rs, 
there arc no helium nuclei and the photodisintc- 
gration term is negligible, while at 7.5 Rs it has 
a value of about Qphoto ~ 10 38 erg/s/cm 2 with 
rapid fluctuations. These fluctuations induce the 
local accretion rate flickering (cf. Fig. [T6|) . on 
a timescale and amplitude much smaller than for 
the viscous instability. 



4.2. Comparison with previous work 

The neutrino dominated accretion flow has al- 
ready been studied in a number of papers, includ- 
ing both 1-D models and multi-D simulations. The 
steady-state 1-D models (e.g Popham et al. 1999; 
Kohri & Mineshige 2002) assumed the disk opti- 
cally thin to neutrinos, and neglected photodis- 
integration cooling. Di Matteo et al. (2002) took 
these two effects into account, and showed that the 
trapped neutrinos dominate the pressure in the 
inner region of the hyperaccreting disc, however 
their equation of state did not include the numeri- 
cal calculation of chemical equilibrium and did not 
incorporate the opacities directly in the EOS itera- 
tions (see also the time-dependent model of Janiuk 
et al. 2004). Kohri, Narayan & Piran (2005) con- 
sidered the neutrino opaque disk and the equilib- 
rium between neutrons and protons and calculated 
the number densities of species by numerically in- 
tegrating their distribution functions. However, 
these authors calculated the gas pressure from the 
ideal gas approximation, and neglected the contri- 
bution of helium to the pressure. In all of these 
papers the disc occurred to be stable against any 
kind of instability. 

On the other hand, in their recent work, Chen & 
Beloborodov (2006) find that the outskirts of the 
disk are gravitationally unstable. The approach 
used by these authors provides a detailed treat- 
ment of the microphysics which is very similar to 
ours; however, some differences between our work 
and theirs must be crucial to the development of 
viscous and thermal instabilities. One difference 
deals with the approximation made for the treat- 
ment of transition region between the neutrino- 
opaque and the transparent matter. In our work, 
we adopted a gray body model, i.e., we introduced 
the b factor to describe the distribution function 
(cf. Eq. 10). This assumption is consistent with 
the two fluids approximation we have made, which 
has recently been studied numerically by Sawyer 
(2003), and shown to be appropriate for the con- 
ditions of these disks. On the other hand, Chen 
& Beloborodov (2006) smoothly connect the opti- 
cally thin and thick regimes by means of interpo- 
lation. A further difference lies in the description 
of the mass fraction of free nucleons. In this work 
we use an expression for X nuc developed by Qian 
et al (1996), while in Chen & Beloborodov (2006) 
X nuc is a function of Y e , which couples the nucle- 
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osynthesis to the electron fraction. 

In our calculations we reach the range of densi- 
ties and temperatures where the nucleons start to 
become partially degenerate. This is accompanied 
by the neutrinos being more and more trapped in 
the gas and helium being destroyed by photodisso- 
ciation. As a result of these calculations, we found 
an additional, unstable branch of solutions for the 
disc thermal balance. 

This supports the recent results of 2-D simu- 
lations by Lee, Ramirez- Ruiz & Page (2005), who 
found the disc opaque to neutrinos to be thermally 
unstable. Their simulations showed that large cir- 
culations develop in the accretion flow. Setiawan, 
Ruffert and Janka (2005) found small fluctuations 
of the accretion rate and neutrino luminosity on 
the dynamical timescale, after the 10-20 msec of 
relaxation period (note, that in our calculations we 
start from the steady-state disc model at a given 
accretion rate, thus having no need for a relaxation 
to the quasi-steady configuration). The equation 
of state used in their work (see also Janka et al. 
1999) is based on the work of Lattimer & Swesty 
(1991). Given the electron fraction, this EOS as- 
sumes the condition of nuclear statistical equilib- 
rium without neutrino trapping, but the evolution 
of the electron fraction is affected by the asym- 
metric neutrino emission from the hot and dense 
matter, which is called 'neutrino leakage scheme'. 
The neutrino leakage scheme focuses on the ef- 
fects of the neutrino trapping on the net neutrino 
emissivities, not on the nuclear statistical equilib- 
rium. The equation of state used in the work of 
Rosswog et al. (2004) is temperature and compo- 
sition dependent, based on the relativistic mean 
field theory (Shcn et al. 1998a, b), and the neu- 
trino cooling is accounted for by the multiflavor 
scheme (Rosswog & Liebendoerfer 2003). 

In our work, we use an equation of state based 
on the (3 equilibrium, including the contribution 
from the trapped neutrino, and neutrino trapping 
effects arc accounted for by the appropriate opac- 
ities. It should be emphasized that most previous 
multi-D simulations neglected the effects of neu- 
trino trapping on the (3 equilibrium, as well as the 
contribution of the trapped neutrinos to the ther- 
modynamical properties of the dense matter. An- 
other difference between our treatment of the EOS 
and the previous numerical simulations is that we 
include the cooling of the photodisintegration of 



helium. Even though the original EOS of Lattimer 
& Swesty (1991) can provide detailed information 
about the composition of the dense matter, this 
information was not considered in order to keep 
the table of the EOS as small as possible (see e.g. 
Ruffert et al. 1996) just for numerical reasons. In 
this way, the disintegration cooling had not been 
investigated without the information on the com- 
position. Our results indicate that photodisinte- 
gration significantly affects the energy balance. 

4.3. Limitations of our model 

We find the thermal- viscous instability to be an 
intrinsic property of the disc for extremely large 
torus densities (about 10 12 g cm" 3 ) and high ac- 
cretion rates (M > 10M Q /s). This is seen both 
in the steady-state results (radial profiles of den- 
sity and temperature) and in the subsequent time 
evolution. 

Thermal and viscous instabilities have been 
studied in the case of standard accretion discs 
around compact objects (Lightman & Eardley 
1974; Pringlc 1977; Shakura & Sunyaev 1976). 
Two main physical processes that lead to disc 
instabilities were invoked to explain the time- 
dependent behavior of various objects: partial 
ionization of hydrogen in the discs of Dwarf No- 
vae (e.g. Meyer & Meyer-Hofmeister 1981; Smak 
1984) and domination of radiation pressure in the 
X-ray transients (e.g Taam & Lin 1984). Such in- 
stabilities do not have to lead to a total disc break- 
down, but rather to a limit-cycle behavior, if only 
an additional (i.e. upper) stable branch of solu- 
tions can be found. This might be a hot state with 
a temperature above 10 4 K, or a slim disc, domi- 
nated by advection (Abramowicz et al. 1988). In 
our 1-D calculation the disc in the GRB central en- 
gine is not stabilized but rather breaks down into 
rings, as no stable solutions are reached (possibly, 
for even higher accretion rates again a stable part 
near the black hole could be formed - but these 
extremely high accretion rates would not be pro- 
duced by any compact merger scenario). There- 
fore, instead of a limit-cycle activity, what we 
find here are several dramatic accretion episodes 
on the viscous timescale. The remaining parts of 
the torus will subsequently accrete and, while ap- 
proaching the central black hole, will get hotter 
and denser, breaking at ~ 7Rs- 

Of course, it would be interesting to study 
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whether such a violent instability would occur also 
in the 2D or 3D simulations. This is indeed likely 
to be the case, since as the multi-dimensional sim- 
ulations of accretion discs show, the instabilities 
derived first in ID are still present in the hydro- 
dynamical simulations of flows with non-Keplerian 
velocity fields (e.g. Agol et al. 2001; Turner 
2004; Ohsuga 2006). Possibly, the instability re- 
gion would be located at other (larger) radii if 
the calculations included the vertical structure of 
the disc: this is dependent on temperature and 
density, which above the meridional plane may be 
larger than the mean value considered in the ver- 
tically averaged model. 

We need to note that our ID calculations do 
not take into account the possible effects of non- 
radial velocity components in the fluid. For exam- 
ple, the inverse composition gradient that leads 
to the disk instability, might be stabilized by ro- 
tation (e.g. Begelman & Meier 1982; Quataert & 
Gruzinov 2000). In the 2-D simulation of Lee et al 
(2005) the neutrino opaque disk exhibits circula- 
tions in the r-z direction. Such meridional circula- 
tions are known to be present in the Kcplerian ac- 
cretion disks (e.g. Siemiginowska 1988), however 
it is unclear if they could always provide a stabiliz- 
ing mechanism for the thermal-viscous instability. 
Possibly, if the nonradial motions of the flow pro- 
vided a stronger stabilizing effect, the disk would 
exhibit oscillations in the viscous timescale, with- 
out breaking, similarly to the outbursts of Dwarf 
Nova disks. 

The assumption of the [3 equilibrium (justified, 
as the mixture of protons, electrons, neutrons and 
positrons is able to achieve the equilibrium con- 
ditions) might also have an effect on this result, 
as the equilibrium conditions reduce the heating 
and entropy in the gas. In fact, the j3 equilibrium 
condition which is satisfied in the innermost part 
of a hyperaccreting disc that is optically thick to 
neutrinos, is /i n = /i p + /i c . Once the disc becomes 
transparent in its outer part, this condition is no 
longer valid. Analytically, it has been derived by 
Yuan (2005) that the condition for (3 equilibrium 
in this case is /i n = fi p + 2/i c . 

4.4. Observational consequences 



lead to a variable energy output on small (millisec- 
ond) timescales. The consequence of this may be 
variability in the gamma ray luminosity, although 
the changes in the local accretion rate may be 
spread by viscous effects (in the lightcurve L v {t) 
integrated over the whole surface of the disc, the 
millisecond variability is somewhat smeared, and 
the amplitudes are not very large). Therefore the 
mass accreted by the black hole may not be vary- 
ing substantially, while some irregularity in the 
overall outflow could help produce internal shocks. 

The thermal-viscous instability, if accompa- 
nied by the disc breaking, may lead to the sev- 
eral episodic accretion events and several re- 
brightenings of the central engine on longer 
timescales, possibly detected in the later stages 
of the evolution. A similar kind of a long-term 
activity is possible also if the disk was not com- 
pletely broken, but exhibited some large accretion 
rate fluctuations on the viscous timescale. 
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Our findings might be relevant for interpret- 
ing some recent observations. The flickering due 
to the photodisintegration of alpha particles may 



19 



A. Appendix 

The neutrino absorption and production rates in the beta processes for all participating particles at 
arbitrary degeneracy have been obtained in the previous works (Reddy, Prakash & Lattimer 1998; Yuan 
2005). In the subnuclear dense matter with high temperatures, the nucleons are generally nondegenerate, 
therefore, the transition reaction rates from neutrons to protons and from protons to neutrons can be 
simplified as follows: 

r p+e -^„+„ c - i' M | 2 / dE e E ePe (E e -Q) 2 f e (l-b e UJ, (Al) 
1 f 00 

r p+ ™+, c = ^\ M \ 2 J dE e E ePe (E e -Q) 2 (l-f e )b e .U c , (A2) 
1 f°° 

r„ +e+ ^ P+Pe = ^\ M \ 2 J dE e E ePe (E e + Qff e+ (l-b e U e ), (A3) 
1 f°° 

r„ +e+ ^ p+Pc - 2^ |M| V dE eEePe(E e +Q) 2 (l-f e+ )b e ,U^ (A4) 

i r Q 

T n ^ p+e - +De = ^\ M \ J dE e E e p e {Q-E e ) 2 {l-f e )(l~b e U c ), (A5) 



r„^ P+e -+ Pe = ^\ M \ 2 J dE e E ePe (Q-E e j 2 f e b e f^. (A6) 

Here Q = (m„ — m p )c 2 , \M\ 2 is the averaged transition rate which depends on the initial and final states 
of all participating particles, for nonrelativistic noninteracting nucleons, \M\ 2 = G|cos 2 #c(l + 3.g A ), here 
Gf — 1.436 x 10 -49 erg cm 3 is the Fermi weak interaction constant, 9c (sin 9q = 0.231) is the Cabibbo angle, 
and <7a = 1-26 is the axial-vector coupling constant. / e ,„ e are the distribution functions for electrons and 
neutrinos, respectively. The "chemical potential" of neutrinos is generally assumed to be zero. The factor 
b e reflects the percentage of the partially trapped neutrinos. When neutrinos completely trapped, b e = 1. 
The corresponding neutrino emissivities for the URCA reactions are given by: 

1 f°° 

q P+ e-^n+, c = ^s^ 1 ! 2 ) dE e E ePe (E e -Qff e (l-b e UJ, (A7) 
1 f°° 

qn+e+^p+ve = ^s\ M \ 2 J dEeE^Ee+Qffe+il-bJ.J, (A8) 

q n ^ P+ e-+p e = Al M ! 2 [ Q dE e E ePe (Q-E e ) 3 (l-f e )(l-b e UJ. (A9) 

The emissivities due to the electron-positron pair annihilation, following the notation of Yakovlev et al 
(2001), is written as: 

« e -+e+^+e 4 - %~ {G+ Vi [8(*lE/2 + ®2Ul) - 2($_ 1 £/ 2 + $2^-l) + 7($o!7i + $iU Q ) 

307T 1 

+ 5(*oy-i + *-i^o)] + 9C^ Wi [*o(^i + U'-i) + (*-i + *i)^o]}, (A10) 

where 

Q c = 3. (^) 9 = 1.023 x 10 23 erg cm- 3 s" 1 , (All) 

C+ Vi = Cy. + C\. and C_„ 4 = C v . — C\., here Cy i and C Ai are the vector and axial- vector constants for 
neutrinos (C Ve = 2sin 2 6l c + 0.5, C Ae = 0.5, Gy M = C Vt = 2sin 2 6> c - 0.5 and C All = C Ar = -0.5). The 
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dimcnsionless functions Uk and $k (k= — 1, 0, 1, 2) in the above equation can be expressed in terms of the 
Fermi-Dirac functions: 

U-i = ^0 2 F 1/2 (r, e ,f3 e ) (A12) 
/2 

Uo = ^/3 3/2 [F 1/2 (r, e ,/3 e ) + /3 e F 3/2 (ry e ,/5 e )] (A13) 
Ui = ^J/3 3/2 [F 1/2 (r /e ,/3 e ) + 2/3 e F 3/2 (r /e ,/3 e ) + /3 e 2 f 5/2 (7 ?e ,/3 e )] (A14) 
U2 = ^/3 3/2 [F 1/2 (r, e ,/3 e ) + 3/3 e F 3/2 (r /e ,/3 e ) + 3/3 2 F 5/2 (r /e ,/3 e )+ / 33F 7/2 (7 ?e ,/5 e )]. (A15) 
Replacing rj e with fy e + in Uk, we get the corresponding expressions for 
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